Jamming transition in granular media: A mean field approximation and numerical 

simulations 
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In order to study analytically the nature of the jamming transition in granular material, we have 
considered a cavity method mean field theory, in the framework of a statistical mechanics approach, 
based on Edwards' original idea. For simplicity we have applied the theory to a lattice model and a 
transition with exactly the same nature of the glass transition in mean field models for usual glass 
formers is found. The model is also simulated in three dimensions under tap dynamics and a jamming 
transition with glassy features is observed. In particular two step decays appear in the relaxation 
functions and dynamic heterogeneities resembling ones usually observed in glassy systems. These 
results confirm early speculations about the connection between the jamming transition in granular 
media and the glass transition in usual glass formers, giving moreover a precise interpretation of its 
nature. 

PACS numbers: 45.70.Cc,05.50.+q,64.70.Pf 



I. INTRODUCTION 

A deep connection between glass transition in molec- 
ular glass formers, structural arrest in colloidal systems, 
and jamming transition in granular media 
has often been stressed in the past few years. In spite of 
the fact that these systems are very different one from 
each other, varying suitably the control parameters, a 
slowdown and a subsequent structural arrest in a solid- 
like disordered state are found in each of them. In [2,01 a 
possible phase diagram for jamming is suggested, which 
takes into account the fact that jamming is obtained ei- 
ther raising the volume fraction or lowering the tempera- 
ture or lowering the applied stress. Colloidal suspensions 
and molecular glass formers are both thermal systems, 
and it is commonly accepted that both colloidal glass 
transition and molecular glass transition are of the same 
type despite of the fact that different control parameters 
may drive the transition. The case of granular materi- 
als is instead very different: They are athermal systems, 
since the thermal fluctuations are significantly less than 
the gravitational energy and the system cannot explore 
the phase space without any external driving. Never- 
theless an exceedingly slowing down is observed when a 
granular material is shaken at low shaking amplitude, 
or flows under a low shear stress, with strong analogies 
with the slowing down observed in^ass formers. Exper- 
imental and numerical studies 0, Q have confirmed 
this connection, however its precise nature is still unclear 

m 

In the present paper in order to study this connection 
we apply a statistical mechanics approach to granular 
media. This approach, which has been extensively devel- 
oped in previous works 0, 0] , is based on an elaboration 
of the original ideas suggested by Edwards . The ba- 
sic assumption is that for a granular system subject to an 
external drive (e. g. tapping), after having reached sta- 



tionarity, time averages coincide with suitable ensemble 
averages over the "mechanically stable" states. We have 
shown that this assumption works for different lattice 
models namely that a generalized Gibbs distribution of 
the stable states describes with good approximation the 
stationary state attained by the system under tapping 
dynamics. Here each tap consists in raising the bath 
temperature to a finite value (called tap amplitude) and, 
after a lapse of time (called tap duration) quenching the 
bath temperature back to zero. By cyclically repeating 
the process the system explores the space of the mechan- 
ically stable states. 

We thus consider one of the above lattice model for 
which the statistical mechanics approach works. The 
model is made up of hard spheres under gravity. Then 
we apply standard statistical mechanics methods in order 
to investigate analytically the existence and the nature 
of a possible jamming transition. More precisely we con- 
sider the Bethe-Peierls approximation using the cavity 
method [ill IT^ : By changing the control parameter a 
phase transition from a fluid to a crystal is found, and, 
when crystallization is avoided, a glassy phase appears. 
The nature of this glassy phase is analogous to that found 
in mean field models for glass formers [T^ . IT^ H^ : In 
particular wc observe a dynamical transition, where an 
exponentially high number of metastable states appears, 
and at a lower temperature a thermodynamic discontin- 
uous phase transition to a glassy state. A brief account 
of these calculations was given in a previous Letter [l5l |. 
We also studied 0| the model in 3d by means of numer- 
ical simulations, and we found that the model under taps 
has a transition from a fluid to a crystal, in a very good 
agreement with he mean field approximation. However 
the numerical simulation was not suitable to study the 
glass transition since the model showed a strong tendency 
towards crystallization. 

For this reason we study here a variant of the model 
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[iS] which has the virtue of avoiding crystaUization. We 
find that the system under gravity evolved by Monte 
Carlo taps presents features characteristic of real granu- 
lar media |l6lll7| . and at low tap amplitudes a dynamical 
transition with properties recalling those of usual glass 
formers. In particular we observe a dynamical non linear 
susceptibility with a maximum at increasing time: This 
behavior, typical of glass formers, is usually interpreted 
as the sign of dynamic heterogeneities in the system. 

In conclusions the results confirm early speculations 
about the deep connection between the jamming transi- 
tion in granular media and the glass transition in usual 
glass formers, giving moreover a precise interpretation to 
its nature. 

In Sect. |n]the mean field phase diagram is discussed. 
The details of calculations are presented in App.slXland 
IbI In particular in App. ^the self-consistency equations 
obtained using the cavity method arc shown. In Sect. IIIII 
the 3d model is presented and the numerical results are 
shown. 



II. MEAN FIELD SOLUTION IN THE 
BETHE-PEIERLS APPROXIMATION 

The model is a monodisperse hard sphere system (with 
diameter V2ao) under gravity, constrained to move on 
the sites of a cubic lattice of spacing gq — 1. The Hamil- 
tonian is given by: 

H^Hhc +rng'^niZi (1) 

i 

where Zi is the height of site i, g = 1 is the gravity 
acceleration, m = I the grain mass, rii € {0, 1} is the 
occupancy variable (absence or presence of a grain on site 
i) and 7i_f/c({"-i}) is the hard core term preventing two 
nearest neighbor sites being simultaneously occupied. 

We have shown in previous papers j^] that the model, 
Eq. J^lj evolving by means of a tap dynamics can be 
described in good approximation by a generalized Gibbs 
distribution of the "mechanically stable" states (i.e. the 
states where the system is found at rest). In particular 
the weight of a given state, {n^}, is: 

e-/3^^({"J).n({nJ), (2) 

where T^onf = I3~^ is a thermodynamic parameter, called 
"configurational temperature" , characterizing the distri- 
bution. The operator n({ni}) selects mechanically sta- 
ble states: n({ni}) = 1 if {n^} is "stable", or else 
n({ni}) = 0. The system partition function is thus the 
foUowing 0: 

2 = ^e-''^«-».n(K}) (3) 

where the sum runs over all microstates, {ni\. 

In the present section we show the phase diagram of 
the model, Eq. (Q, obtained using a mean field theory 



in the Bethe-Peierls approximation (see a-nd ref.s 

therein), based on a random graph (plotted in Fig. ^ 
which keeps into account that the gravity breaks up the 
symmetry along the z axis. This lattice is made up by 
H horizontal layers (i.e., z S {l,...,i?}). Each layer is 
a random graph of connectivity, fc — 1 = 3. Each site in 
layer z is also connected to its homologous site in z — 1 
and z + 1 (the total connectivity is thus fc -|- 1). Locally 
the graph has a tree-like structure but there are loops 
whose length is of order In TV, insuring geometric frustra- 
tion. In the thermodynamic limit only very long loops 
are present. The details of calculations are given in ap- 
pendiceslXlandlBl fsee also Ref.s where this mean 

field theory was first introduced). 




z+1 



z 



z-1 



FIG. 1: In the mean field approximation, tlie grains are lo- 
cated on a Bethe lattice, sketched in the figure, where each 
horizontal layer is a random graph of given connectivity. Ho- 
mologous sites on neighboring layers are also linked and the 
overall connectivity, c, of the vertices is c = -|- 1 = 5. 

We solve the recurrence equations found in the Bethe- 
Peierls approximation in three cases: 1) A fluid- like ho- 
mogeneous phase; 2) a crystalline-like phase character- 
ized by the breakdown of the horizontal translational in- 
variance; 3) a glassy phase described by a 1-step Replica 
Symmetry Breaking (IRSB). The details of the calcula- 
tions are shown in Appendices. 

The results of the calculations are summarized in 
Fig. 121 where the bulk density at equilibrium, <f> = 
Ns/{2{z) — 1) ^3 (where (2) is the average height) is 
plotted as a function of the configurational temperature, 
Tconf, for a given value of the number of grains per unit 
surface, Ng. We found that at high Tconf a homogeneous 
solution corresponding to the fluid-like phase is found. 
By lowering Tconf at T^ a phase transition to a crystal 
phase (an anti-ferromagnetic solution with a breakdown 
of the translation invariance) occurs. The fluid phase still 
exist below Tm as a metastable phase corresponding to 
a supercooled fluid when crystallization is avoided. Fi- 
nally a IRSB solution (found with the cavity method 
characterized by the presence of a large number of 
local minima in the free energy appears at To, and 
becomes stable at a lower point Tk, where a thermody- 
namic transition from the supercooled fluid to a IRSB 
glassy phase takes place. The temperature Td, which is 
interpreted in mean held as the location of a dynamical 
transition where the relaxation time diverges, in real sys- 
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tems might instead correspond to a crossover in the dy- 
namics (see [Hill Hi and Ref.s therein). <^{T^onf) has 
a shape very similar to that observed in the "reversible 
regime" of tap experiments [Til I2H . The location of the 
glass transition, Tk, corresponds to a cusp in the func- 
tion ^{Tconf)- The dynamical crossover point Td might 
correspond to the position of a characteristic shaking am- 
plitude r* found in experiments and simulations where 
the "irreversible" and "reversible" regimes approximately 
meet. 
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FIG. 2: The density, <1> = N,/{2{z) - 1), for A^^ = 0.6 as a 
function of Tconf- ^max is the maximum density reached by 
the system in the crystal phase. 

In Fig.|31the phase diagram obtained by varying Ng is 
shown. The dashed vertical line in figure corresponds to 
the value of Ng chosen in Fig. |2 
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FIG. 3: The system mean field phase diagram is plotted in 
the plane of its two control parameters {Tconf ,Ns). 

The model, Eq. simulated in 3d by means of Monte 



Carlo tap dynamics '15| presents a transition from a fluid 
to a crystal as predicted by the mean field approximation, 
density profiles in good agreement with the mean field 
ones, and in the fiuid phase a large increase of the relax- 
ation time as a function of the inverse tap amplitude. In 
the following section we study a more complex model for 
hard spheres, where an internal degree of freedom allows 
to avoid crystallization 13] . 

III. HARD SPHERES WITH AN INTERNAL 
DEGREE OF FREEDOM 

The Hamiltonian of the model is 

ninj(j>ij{ai,aj) + rng'^^^riiZi, (4) 

{i]) 

where Zi is the height of site i, g — 1 is the gravity 
acceleration, m = 1 the grain mass, rii g {0, 1} is the 
occupancy variable (absence or presence of a grain on 
site i), (Ji e {1, . . . , 9} represents the internal degree of 
freedom (which we call spin), and <j)ij{(Ji, aj) is the inter- 
action energy between spins. Different values of the spin 
correspond to different positions of the particle inside the 
cell. It is reasonable that a few number of internal states 
might be enough to catch the main features of real sys- 
tems. 

As in Ref. study a simple realization of the 

model described by Eq. Q. Interpreting the spin as po- 
sition of the particle in the cell, our choice can be easily 
visualized in 2d, as shown in Fig. ^ We partition the 
space in square cells, and subdivide each cell into four 
internal positions (namely q — 4). When a cell is oc- 
cupied by a particle in any given position, a hard-core 
repulsion excludes the presence of particles in some of 
the internal states of the neighboring cells (namely the 
interaction 0^ [cri , aj ) is chosen zero if the positions ai 
and aj are "compatible", and infinite otherwise). This 
choice can be interpreted as a coarse grained version of 
a hard sphere system in the continuum. In 3c? we subdi- 
vide the space into cubic cells, and considers six internal 
positions instead of four. 
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FIG. 4: The model in two dimensions: the space is partitioned 
in square cells, and each cell can be occupied by at most one 
particle in anyone of the four shown positions (little circles). 
A particle in any given position (large shaded circle) excludes 
the presence of particles in any of the black colored positions. 
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FIG. 5: The bulk density, -I- = N/L^{2{z) - 1), is plotted 
as function of Tr for tq = 10 MC steps /particle. The empty 
circles correspond to stationary states, and the black stars 
to out of stationarity ones. $max is the maximum density 
reached by the system in the crystal phase, $max = 6/7. 
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FIG. 6: The density profile, cr(z), as function of the height, 
z, for Tr = 0.20 and tq = 10 M C steps / particle. 



In the Monte Carlo simulations, N = 433 grains are 
confined in a 3c? box of linear size L = 12 (i.e. Ng — 3), 
between hard walls in the vertical direction and with pe- 
riodic boundary conditions in the horizontal directions. 
We perform a standard Metropolis algorithm on the sys- 
tem. The particles, initially prepared in a random con- 
figuration, are subject to taps, each one followed by a re- 
laxation process. During a tap, for a time tq (called tap 
duration), the temperature is set to the value Tp (called 
tap amplitude), so that particles have a finite probabil- 
ity, Pup ~ e"™^/"^"", to move upwards. During the re- 
laxation the temperature is set to zero, so that particles 
can only reduce the energy, and therefore can move only 



downwards. The relaxation stops when the system has 
reached a blocked state, where no grain can move down- 
wards. Our measurements are performed at this stage 
when the shake is off and the system is at rest. The 
time, t, is the number of taps applied to the system. 

In the following the tap duration is fixed, tq = 
lOMC steps /particle, and different tap amplitudes, Tr, 
are considered. In Fig. [S] the bulk density, $ = 
N/L^{2{z) - 1), is plotted as a function of Tp: $(Tr) 
has a shape resembling that found in the "reversible 
regime" of tap experiments |l6l I2H , and moreover very 
similar to that obtained in the mean field calculations 
and shown in Fig. [5] At low shaking amplitudes (cor- 
responding to high bulk densities) a strong growth of 
the equilibration time (i.e. the time necessary to reach 
stationarity) is observed, and for the lowest values here 
considered (the black stars in Fig.O the system remains 
out of stationarity. In this region the density profile, 
a{z) = ni{z) (where the sum is done over the 

sites i in the layer z = 1, . . . , L H^), is almost constant 
until a given layer and sharply decays to zero (see Fig. 
EJ, as found in real granular media p?7| . In conclusions 
the system here studied presents a jamming transition at 
low tap amplitudes as found in real granular media. 

In order to test the predictions of the mean field cal- 
culations, in the following we measure quantities usually 
important in the study of glass transition: The relaxation 
functions, the relaxation time and the dynamical suscep- 
tibility, connected to the presence a dynamical correla- 
tion length. 

In particular we calculate the two-time autocorrelation 
functions: 

C{t,t^) = ■^'^ni{t)ni{ti^)ai{t)-ai{ty,), (5) 

i 

where ai are unit length vectors, pointing in one of the 
six coordinate directions, representing the position of the 
particles inside the cell; the average (...) is done over 
16 — 32 different realizations of the model obtained vary- 
ing the random number generator in the simulations, and 
the errors are calculated as the fluctuations over this sta- 
tistical ensemble. For values of t^^ long enough, the sys- 
tem reaches a stationary state, where the time translation 
invariance is recovered, i.e., C{t, ty^) — C{t — t^). In this 
time region, by averaging C(t' ,tyj) over t' and such 
that t = t' — is fixed, we calculate the "equilibrium" 
autocorrelation functions 



{q{t)) = {C(t'~t^)), 
and the dynamical non linear susceptibility 
X{t) = {q{i))\ 



(6) 



(7) 



As shown in Fig. [7| at low values of the tap amplitudes, 
Tr, two-step decays appear, well fitted in the interme- 
diate time region, by the /3— correlator predicted by the 
mode coupling theory for supercooled liquids pMl24j (the 



5 




FIG. 7: The "equilibrium" autocorrelation function, 
(g(t)), plotted as function of t, for tap amplitudes 
Tr = 0.60, 0.50, 0.425, 0.40, 0.385, 0.365, 0.36 (from bot- 
tom to top). The continuous line in figure is the /3— correlator 
of the mode coupling theory with exponent parameters 
a = 0.30 and b — 0.52. The dashed line is a stretched 
exponential oc exp[~{t/r)^] with f3 = 0.70. 

continuous curve in Fig. [7)), and at long time by stretched 
exponentials (the dashed curve in figure). The relaxation 
time, T, is defined as (q'(r)) ~ 0.1. 

In Fig.|Slthe relaxation time, r, is plotted as a function 
of the density, $. As found in many glass forming liquids, 
t($) is well fitted by a Vogel-Fulcher for the entire range, 
even if we can identify a first region where t($) is fitted 
with good approximation by a power law. The power 
law divergence can be interpreted as a mean field behav- 
ior, followed by a hopping regime. Note that the model, 
Eq. Q), studied in absence of gravity by means the usual 
Monte Carlo Metropolis , exhibits a divergence of the 
relaxation time as a power law, and no crossover to a hop- 
ping regime is observed. We suggest that in the present 
case the tap dynamics favors the equilibration via hop- 
ping precesses. 

In Fig.l^lthe relaxation time, t, is plotted as a function 
of the tap amplitude, Tr: A clear crossover from a power 
law to a different regime is again observed around a tap 
amplitude To, corresponding to the value of the density, 
<I>(Td) ~ $£1, where a similar crossover has been found 
in Fig. El 

The divergence of the relaxation time at vanishing tap 
amplitude is consistent with the experimental data of 
Philippe and Bideau |3 and D'Anna et al. 4]. Their 
findings are in fact consistent with an Arrhenius behav- 
ior as function of the experimental tap amplitude inten- 
sity. However a direct comparison with our data is not 
possible since we do not know the relation between the 
experimental tap amplitude and the tap amplitude in our 
simulations. A more direct comparison would be possible 
if the experimental data were plotted as function of the 
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FIG. 8: The relaxation time, r, as function of the bulk density. 
The continuous line is a Vogel-Fulcher, e^'''*''"*' ,with 
$c = 0.81 ± 0.01 and A = 0.49 ± 0.10. The dashed line 
is a power law, ($d — 'I>)~^S with <I>i3 — 0.76 ± 0.01 and 
71 = 2.04 ±0.10. 
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FIG. 9: The relaxation time, r, as function of the tap 
amplitude inverse, Tj7^. The dashed line is a power law, 
(Tr - Td) with Td = 0.40±0.01 and 72 = 1.52±0.10. The 
continuous line is an Arrhenius fit, e^^^r^ .j^j^h A = 17.4±0.5 
(the data in this region are also well fitted by both a super- 
Arrhenius and Vogel-Fulcher laws). 



bulk density, as we did in Fig. |S| 

The dynamical non linear susceptibility, x{t), plotted 
in Fig. llOl at different Tr, exhibits a maximum at a time, 
t* (Tr ) . The presence of a maximum in the dynamical non 
linear susceptibility is typical of glassy systems . 
In particular the value of the maximum, x(i*), diverges 
in the p-spin model j2||| as the dynamical transition is ap- 
proached from above, signaling the presence of a diverg- 
ing dynamical correlation length. In the present case the 
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value of the maximum increases as Tp decreases (except 
at very low Tp where the maximum seems to decrease 
[231 )■ The growth of in our model suggests the 

presence of a growing dynamical length also in granular 
media. 
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FIG. 11: Three kinds of branches exist here: 1) "Side" branch: 
the root site is connected to fc — 2 neighbors on its layer, one in the 
upper and one in the lower layer; 2) "Up" branch: the root site is 
connected to fc — 1 neighbors on its layer and one in the upper layer; 
3) "Down" branch: the root site is connected to fc — 1 neighbors on 
its layer, one in the lower layer. 
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FIG. 10: The dynamical non linear susceptibility, (nor- 
malized by x(^o), the value at to = 1) as a function of t, for tap 
amplitudes Tr = 0.60, 0.50, 0.425, 0.41, 0.40, 0.385, 0.3825 
(from left to right). 



IV. CONCLUSIONS 

In conclusions using standard methods of statistical 
mechanics we have investigated the jamming transition in 
a model for granular media. We have shown a deep con- 
nection between the jamming transition in granular me- 
dia and the glass transition in usual glass formers. As in 
usual glass formers the mean field calculations obtained 
using a statistical mechanics approach to granular media 
predict a dynamical transition at a finite temperature, 
Td, and, at a lower temperature, Tk, a thermodynamics 
discontinuous phase transition to a glass phase. In finite 
dimensions 1) the dynamical transition becomes only a 
dynamical crossover as also found in usual glass formers 
p2. 0, ll^l (here the relaxation time, r, as a function 
of both the density and the tap amplitude, presents a 
crossover from a power law to a different regime); and 2) 
the thermodynamics transition temperature, Tk, seems 
to go to zero (the relaxation time, r, seems to diverge 
only at Tp ~ 0, even if a very low value of the transition 
temperature is consistent with the data). 
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APPENDIX A: MEAN FIELD SOLUTION 

We consider the Hamiltonian, Eq. plus a chem- 
ical potential term which controls the overall density. 
Hard Core repulsion prevents two connected sites to be 
occupied at the same time. We adopt a simple defini- 
tion of "mechanical stability": a grain is "stable" if it 
has a grain underneath. For a given grain configura- 
tion r = {ui}, the operator 11^ has a simple expres- 
sion: Ur = limK^oc e^P {-KTiEdw}, where TiEdw = 
J2iSni{z),iSniiz-i),oSni(z-2),o (for clarity, we have shown 
the z dependence in ni{z)). 

The random graph, Fig.^ keeps into account that the 
gravity breaks up the symmetry along the z axis. This 
lattice is made up by H horizontal layers |2^ occupied 
by hard spheres (two numbers identify a site of the lat- 
tice: The height of the layer, z G {1,...,H}, and the 
position in the layer, i). Each layer is a random graph 
of connectivity, fc — 1 = 3. Each site in layer at height 
z is also connected to its homologous site in z — 1 and 
z + 1 (the total connectivity is thus fc + 1). The local 
tree-like structure of the lattice allows to write down it- 
erative equations a la Bethe, where the partition function 
of each site is written in terms of the partition functions 
of the neighbor sites. We have to introduce the concept 
of "branch": a branch is a graph where a root site, i, 
has only k neighbors. In the present case three kinds of 
branches exist (see Fig. ^3 : "up" (resp. "down") branch 
where the root site has fc — 1 neighbors on its same layer 
and one in the upper (resp. lower) layer; and "side" 
branch where the root has fc — 2 neighbors on its layer, 
one in the upper and one in lower layer. 

(i z) (i z) 

Define Zq ^ and Z{ ^ the partition functions of a 
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"side" branch with root site i at height z restricted re- 
spectively to configurations in which the site i is empty 

or filled by a particle. Z^^]^'^ and Zg]^'' (resp. Z^^^) are 
the partition functions of the "up" branch restricted re- 
spectively to configurations in which the site i is filled 
by a particle, or empty with the upper site filled (resp. 

empty). Finally z'^l"^^ and z'^^'^ (resp. Zq '^'') are the 
partition functions of the "down" branch restricted re- 
spectively to configurations in which the site i is filled 
by a particle, or empty with the lower site empty (resp. 
filled). 

The partition function of the branch ending in site i 
can be recursively written in terms of the partition func- 
tions of the neighbor sites. Summing over all the possible 
configurations of the neighbor sites, we obtained that the 
partition function of a "side" branch with root site i at 
height z is: 
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In the same way we can write the recursion relations for 
the "up" branch: 
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Finally for the "down" branch we have: 
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In the following we consider the limit K 00 in order to 
take into account the constraint on the mechanical sta- 
bility. It is convenient to introduce five local "cavity" 
fields on each site /li*'^'', hu'^\ gu'^\ h^d^^ 5d '^^ 
fined by the following relations: e'^'** ' ' — z'"^'^^ / z''^'^^ \ 

r,(i,z) i-TyKi^z) nJ-^-'l r,{i,z) i-7^(i,z) ^ 

/^o,d ; e'^J^'^ = /^o,u ■ In these new vari- 
ables the recursion relations are more easily written: 



C0K 



^p{p-mgz) 



fe-2 

J=l 



^Pin-mgz) 



g/3/i<^-)-j-l 

(i.^-l) 



k-l 

i=i 



(Al) 



{l + e 



P9)l 



ah'; 



(A4) 



^(3(p-mgz) -ph'-^'' 



k-l 

J=l 



g/3/i<^'^>-)-l 



Note that in the case k = 1 the problem reduces to a 
simple one-dimensional chain: In this case the recursive 
method is equivalent to the transfer matrix method and 
gives the exact solution. 

From the iterative solution of Eq.s (|A4p it is possible to 
compute the system free energy. Generalizing the proce- 
dure followed in we calculate the free energy density, 
F, in the thermodynamic limit from the variation of the 
free energy going from a random graph with H layers 
and TV sites on each layer to one with H layers and N + 2 
sites on each layer. In order to do that we define the 
following intermediate object: a random graph with H 
layers and N sites in each plane such that 2(k + 1) sites 
on each plane are connected only to k neighbors instead 
of fc -|- 1 . In particular on each layer 2 sites are not con- 
nected to sites on the higher layer ("down" branches), 2 
sites are not connected to sites on the lower layer ("up" 
branches) and the other 2{k — 1) are connected only with 
k — 2 sites in the plane instead of fc — 1 ("side" branches). 
From this intermediate object a random graph with H 
layers and + 2 sites on each layer (all connected to 
k + 1 sites) can be obtained adding 2 new sites to each 
plane and connecting each of the new sites with fc — 1 
side branches on their respective planes, one up branch 
in the upper layer and one down branch in the lower layer 
(see Fig. I12|l . This operation is called "site addition" . A 
random graph with H layers and N sites on each layer 
(all connected to fc -I- 1 sites) is instead obtained from the 
intermediate object adding for each layer 2 links between 
the up branches at height z and the down branches at 
height z — 1, and (fc — 1) links between the sides branches 
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on each layer (see Fig. ll2|) . This operation, which allows 
to saturate all the missing links, is called "link addition" . 

Therefore the variation of the free energy when going 
from NH to {N + 2)H sites (i.e. a random graph with 
two sites more on each layer) is related to the free energy 
shifts (see Fig.^J for a site addition (AF^^^^) and for two 
different kinds of link addition (Af/^^\ ^ and Af/^^ ^- 



link. 2 ' 



H H 

F^+,-F^ = 2^A^^it-(fc-l)E^^/42 

2 = 1 Z=l 

— 9 A p(^) 



Site 

Addition 




Link 
Addition (2) 



Link 
Addition (1) 

O z+1 



6 z 
d 



s O 



-O s 



where Fn+2 ~ Fn is obtained as (i^7v+2 — Po) — (Pn — 
Fo), and Fq is the free energy of the intermediate object 
described above. We assume that in the thermodynamic 
limit the free energy is linear in A'^. The free energy 
density is then: 



" {z} (fc- 1) V- A 

site c 



H-l 



x^af'-^' -^af'-'-' 



7. = \ 



z=l 



(A5) 

In terms of the local fields the free energy shifts due to 
the addition of a site i at height z reads: 



fc-i 

n(' 



"i 



(A6) 

') 



The free energy shift due to a link addition between a 
down branch with the root site at height z and an up 
branch with the root site at height z + 1 is given by: 



+ e 



(A7) 



Finally the free energy shift due to a link addition be- 
tween two side branches with root sites i and j at height 
z is: 



1 



+ e 



f3h' 



(A8) 



In order to compute the free energy of the system we have 
to compute the mean values of the free energy shifts for 
all the sites at a given height and for all the possible real- 
ization of the lattice. In the following these mean values 
will be computed in three different cases: 1) A fluid-like 
homogeneous phase; 2) A crystalline-like solution charac- 
terized by the breakdown of the horizontal translational 
invariance; 3) A glassy phase by a 1-step Replica Sym- 
metry Breaking. 

The fluid-like solution is obtained by setting that the 
local fields on each layer are the same for all sites of the 



FIG. 12: Site addition: a new central site at height z is connected 
to A: — 1 side branches (s) with the root sites in the same layer, to 
one up (u) branch with the root site in the higher layer and to one 
down (d) branch with the root site in the lower layer; Link addition 
(1): a link between a down branch with the root site at height z 
and an up branch with the root site in the higher layer is added; 
Link addition (2): a link between two side branches with the root 
site in the same layer is added. 



layer ({h^'^^)} = {h(^)} Vi). In this case Eq.s be- 
come 5H — 1 algebraic coupled equations and they are 
easily solved finding the fixed points. This homogeneous 
(Replica Symmetric) solution is characterized by horizon- 
tal translational invariance and is found to be stable for 
high values of the configurational temperature, Tconf, or 
for low values of the number of grains per unit surface, 
Ns- In this case the free energy is easily computed from 
Eq.s (jA6|) . (jA8|l and (|A7|) . since in this case all the quan- 
tities are site independent. From the free energy, F, we 
derive the density profile, a{z) = {ni{z)): 



a{z) 



(A9) 



the number of particles per unity of surface, Ng = 
^^-^ct(z), and the gravitational energy density, E = 
I]f=i rngzaiz). From the relation F = E - TS - fj,Ns, 
we also calculate the entropy per lattice site, S — ~l3F — 
P^iNs + pE. 

In the crystalline (Replica Symmetric) solution the lo- 
cal fields are different on different sites (breakdown of 
translational invariance), but do not fluctuate from site 
to site. This is achieved by the introducing two sub- 
lattices, a and 6, and different local fields on each lattice. 
The merging is done taking into account the structure of 
the crystalline phase. In our case, each site of the sub- 
lattice a (resp. 6) is connected with fc -I- 1 sites of the 
sub-lattice b (resp. a). The crystal periodicity is thus 
two lattice spacings. Schematically Eq.s (|A4p for each 
layer become: 



{ha}^ 

{hb} 



f(/3,M,{hb}) 
f(/3,/i,{ha}) 
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where {ha} and {hb} are the sets of all local fields re- 
spectively on the two sub-lattices. This is a system of 
2{5H — 1) algebraic coupled equations. The free energy 
is computed from the fixed points of these equations. For 
a given Ns, by lowering Tconf, a phase transition from the 
fluid to the crystal is found at (see Fig.Oll. 

The fluid phase still exist below r,„ as a metastable 
phase corresponding to a supercooled fluid when crys- 
tallization is avoided. Nevertheless, the entropy per site 
predicted by the fluid solution becomes negative when 
the temperature is lowered, or the packing fraction is 
increased. Hence, this solution is not appropriate to de- 
scribe the high Ns or low Tconf region. A solution char- 
acterized by the presence of a large number of local min- 
ima of the free energy is found in this region. In this case 
the local flelds may fluctuate. To describe this situation 
we have to introduce three probability distributions on 
each layer 'P"^(/i«, ffu), V- ^Qis) and Vf.,{hd,gd) deflned 
as the probability of finding the flelds hu'^'' and (/u (or 

respectively h^s'^\ or h^^^'^ and g^^^^) on site i at height 
z equal to /i„ and (or respectively to hs, or to hd and 
gd)- Since the glassy phase is expected to be transla- 
tional invariant, we work in the factorized case in which 
the probability distributions at a given height are equal 
for aU the sites of the layer (P";^''' = V^'"^'^). 

Within the one-step Replica Symmetry Breaking 
ansatz of the cavity method (see appendix IB|| the recur- 
sion relations for the flelds are replaced by self consistent 
integral equations for the probability distribution of the 
flelds. For the "up" merging the self consistent integral 
equation reads: 



v:{K,gl) 



k-l 

i=i 



(AlO) 



i,z+l) 



mAF};' 



X S{K-h^:'^-^)Sigi-gi^'^^)e-P 

where Ci is a constant insuring the normalization of 'P", 
hu'^'' and (?« are the local fields defined by Eq.s l|A4|l . 
m € [0, 1] is the usual IRSB parameter to be obtained 
by the maximization of the free energy with respect to 
it, and AF^^ is the free energy shift in the "up" merging 
process. This quantity is computed by using that the 
addition of a site i at a certain height z ( "site addition" ) 
is the result of an "up" merging process, which bring to 
a new "up" branch with root site i, plus a link addition 
between this branch and a down branch at height z — 1: 



(i,z-l|z) 



'-^-'^^ site — '-^^up ^ '-^-'^^Hnk,! 

From this equation we obtain that: 



O.u 



-^0,11 lli=l ^0,s 



(All) 



(A12) 



From Eq.s (|A1IA3|) the free energy shift AFup^' has 
simple expression in terms of the local fields. 



In the same way we can determine the self consistency 
equations for the other two kinds of merging: 



(i,z-l),(i,z-l)^d 



dhr-'dg^'^'VUiK 



(i.z — l) (i.z 

^9d 



(A13) 



and 



VtiK^g^d) = C3 / n [dh^r^Vlihi^^^'^)] (A14) 

X 'd4''-''di-'^VtM''-'^.^''~'')_ 
X 5(/,^-4-))J(g|-5(-))e-^™AF<;-l^ 

For the "side" and the "down" merging one has that: 

= ^F^^ + ^Fl:^-l (A15) 

and 



This yields to: 
and 



i^.z) 



z. 



(i,z) 
0,s 



-7^{i,z+l)—(i,z-l) r^U,z) ' 



0,d 



1=1 



(A16) 



(A17) 



-y{i,Z + \) fc_l (-,-,2) ■ 



(A18) 



For any value of /3, /i and m we solve Eq.s (jSini, (Pl3)l 
and ljA14(l iteratively, discretizing the probability distri- 
butions until the whole procedure converged. 

From the probability distributions we compute the free 
energy density of the system: according to Eq. ljA5|l we 
have to find the average values of the free energy shifts 
due to link and site additions. The free energy shift due 
to site addition is given by: 



-/3mAF,it,(z)\ ^ 



n \d}^^^vi{}^^^'^^) 

X [dh^-'^dg^-'-'^VUih, , 



(A19) 

(i,z-l) (i,z-l) 



X e 



-f3mAFsitcii,z) 
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For the first kind of link addition we have: 



(A2o) 

Finally for the second kind of link addition we find; 



(e-''™^^--> = / n \dh^'''vnh^-^y: 



g link, 2 



In the previous relations AFsUeih z)^ 



{i,z-l\z) 
link.l 



(A21) 



and 



are function of the local fields according to 
Eq. Unj, (HH) and (|X8l) . 

The total free energy density of the system is, accord- 
ing to Eq. (|X5)l : 



Ffml = - 



/3m 



^loge-^"^^"- 



(A22) 



H 



E 

H-1 



ige 



The parameter m is fixed by the maximization of the 
free energy with respect to it. The justification for that 
is in the replica method, since m turns out to be the 
breakpoint in Parisi's order parameter function at the 1- 
step RSB level. For a spin glass it has been rigorously 
proved that in the limit k oo, F[m] is a lower bound 
to the correct free energy, so it is natural to find the 
preferred value of m by the maximization of i^[m]. 



APPENDIX B: SELF-CONSISTENCY 
EQUATIONS IN THE CAVITY METHOD 

In this appendix we show how to obtain the self- 
consistency integral equations HA10|) . (|A13p and (IA14|) 
using the cavity method in the 1-step RSB ansatz 11, 12]. 
The region at high packing fraction (or at low configura- 
tional temperature) is characterized by the existence of 
many pure states. Let J^{F) the number of pure states 
for a given value of the free energy of the system. The 
fmiction S(i^) = log J\f{F) is called complexity. We as- 
sume that within one pure states a the local fields hu'a , 
^d'a ^^'i 9d'a^ different cavity sites are 
uncorrelated. Therefore Eq.s ljA4|l continue to hold in 
any given pure state. In this case we have to make a 
statistical description of the solutions of Eq.s ljA4p in the 



different pure states, taking into account the number of 
pure states for a given value of the free energy. 

Let us consider, for example, the "up" merging of k 
cavity sites in a site i at height z. As said before, in each 
pure state a the local fields in the k cavity sites are not 
correlated. Nevertheless in each pure state a the local 

fields, {hu^a\ gu.'a), and the free energy shift, AFu,a\ 
due to the merging are correlated, since they are both 
functions of the local fields in the neighbor sites in the 
state a, according to Eq.s (|A4I) and (|Alip . Let us define 
Sz{hl^, g^, AF^) as the probability distribution of find- 
ing the fields {hu^a , gu,a^) and the free energy shift AF^ 
after an up merging at height z. Because of the recur- 
sion relations of Eq.s ljA4p and IjAlip . this distribution 
probability has to verify the following iteration relation: 



fc-i 



(Bl) 



s.{h:,g:,AF:) = / n 

X S{h:~h^:-^)6{g:-gi^'^^ 
X 5{AF:-AFi^'^^). 



In order to determine the probability distribution for the 
local fields self-consistently we have to make the integra- 
tion over all possible free energy shifts: 



J d{AF:)S.{h:, g:, AF^WiF - AF) 
d(AF„-)5.(/iS,5^,AFJ)e^(^-^^"). 



Since we are interested only in the local minima with the 
lowest free energies, we expand the exponent to the first 
order in AF^: 

'P:{h:,gl)^C^ J d{AF:)S.{h:,gl,AF:)eM-PmAF:), 



where the parameter m £ [0, 1] is: 

1 9E 



(B2) 



(B3) 



and Ci is a normalization constant. Actually the first 
order expansion means that the density of pure states 
for a given value of the free energy is A/" ~ exp(m(-F — 
Fref)), where Fref is a reference free energy whose value 
is completely irrelevant. This form of the density of states 
is the same found in the 1-step RSB formulation. 

By integrating over AF^, the delta function selects 
only the right value of the free energy shift given in 
Eq. l|IT2|) . We thus have: 



(B4) 
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the same way it is possible to obtain the equations for 
the "side" and the "down" merging. 



We have thus obtained the self consistency Eq. (|A10|) . In 
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